# refit models from last time
plants <- read.table("ecol 563/jimsonweed.txt", header=T, sep='\t')
plants[1:8,]
 
# fixed effects RCBD model
mod2 <- lm(lw.rat~factor(pot)*type, data=plants)
anova(mod2)
mod1 <- lm(lw.rat~factor(pot)+type, data=plants)
anova(mod1)
summary(mod1)
 
# mixed effects RCBD model using nlme
library(nlme)
mod2.lme <- lme(lw.rat~type, random=~1|pot, data=plants)
anova(mod2.lme)
anova(mod1)
summary(mod2.lme)
coef(mod1)
fixef(mod2.lme)
ranef(mod2.lme)
coef(mod2.lme)
predict(mod2.lme)[1:12]
predict(mod2.lme, level=0)[1:12]
 
# mixed effects RCBD model using lme4
library(lme4)
mod2.lmer <- lmer(lw.rat~type+(1|pot), data=plants)
anova(mod2.lmer)
summary(mod2.lmer)
# there is no predict method for lmer objects
predict(mod2.lmer)
class(mod2.lmer)
 
#determine what functions have special predict methods
methods("predict")
 
# fitted works for both lmer and lme models
fitted(mod2.lmer)[1:12]
fitted(mod2.lme)[1:12]
 
# assemble estimates in data frame with raw data
new.dat3 <- data.frame(plants, fix.ests=predict(mod1), mix.ests=fitted(mod2.lmer))
new.dat3[1:8,]
fixef(mod2.lmer)
# add population-average mean to data frame
new.dat3$pop.mean <- ifelse(new.dat3$type=='G', fixef(mod2.lmer)[1], fixef(mod2.lmer)[1]+fixef(mod2.lmer)[2])
new.dat3[1:8,]
 
library(lattice)
# dot plot of raw data
dotplot(factor(pot)~lw.rat|type, data=new.dat3)
# dot plot of raw data using a panel function
dotplot(factor(pot)~lw.rat|type, data=new.dat3, panel=function(x,y) {
panel.grid(v=0, h=-16)
panel.xyplot(x, jitter(as.numeric(y)), col=1, cex=.5, pch=16)})
 
# dot plot of raw data plus treatment and block means
dotplot(factor(pot)~lw.rat|type, data=new.dat3, panel=function(x,y,subscripts) {
panel.grid(v=0, h=-16)
panel.xyplot(x, jitter(as.numeric(y)), col=1, cex=.5, pch=16)
panel.points(new.dat3$fix.ests[subscripts], y, pch=1, col=4)
panel.points(new.dat3$mix.ests[subscripts], y, pch=8, col=2)
panel.abline(v=new.dat3$pop.mean[subscripts], lty=2, col='seagreen')
},scales=list(x='free'))
 
# RCBD with no replication of treatments within blocks
library(faraway)
data(oatvar)
oatvar[1:8,]
dim(oatvar)
table(oatvar$block,oatvar$variety)
# additive model
g <- lm(yield~block+factor(variety), data=oatvar)
anova(g)
# interaction model does not yield standard error estimates
g2 <- lm(yield~block*factor(variety), data=oatvar)
anova(g2)
summary(g2)
 
# Tukey test for non-additivity
coef(g)
blockcoefs <- c(0,coef(g)[2:5])
trtcoefs <- c(0,coef(g)[6:12])
oatvar$alpha.gamma <- rep(trtcoefs, each=5) * rep(blockcoefs,8)
oatvar[1:8,]
g1 <- update(g, .~. +alpha.gamma)
anova(g1)
 
# graphical checks for block x treatment interaction
dotplot(block~yield, group=variety, data=oatvar, pch=1:8)
interaction.plot(oatvar$variety, oatvar$block, oatvar$yield)
 
# complicated ANOVA models with two kinds of blocks
nitro <- read.csv('ecol 563/nitro.csv')
dim(nitro)
nitro[1:8,]
table(nitro$lh, nitro$func, nitro$n, nitro$p)
 
#create a treatment variable
nitro$trt <- paste(nitro$lh, nitro$func, nitro$n, nitro$p, sep='.')
nitro[1:8,]
 
# not all treatments are represented in each block
table(nitro$trt, nitro$block)
table(nitro$trt, nitro$phy)
# response varies by phy
boxplot(pN~phy, data=nitro)
 
# 4-factor ANOVA model with blocks
mod3 <- lm(pN^2 ~ factor(block)+factor(phy) + factor(lh)*factor(n)*factor(func)*factor(p), data=nitro[nitro$tag!=444,])
# 6-factor interaction model fails to estimate
mod0 <- lm(pN^2 ~ factor(block)*factor(phy)*factor(lh)*factor(n)*factor(func)*factor(p), data=nitro[nitro$tag!=444,])
anova(mod0)
 
# compare Type I and Type II tests
anova(mod3)
library(car)
Anova(mod3)
 
# drop 4-factor and 3-factor interactions
mod3a <- lm(pN^2 ~ factor(block)+factor(phy) + (factor(lh)+factor(n)+factor(func)+factor(p))^2, data=nitro[nitro$tag!=444,])
anova(mod3a)
# add back in 3-factor interaction with lowest p-value
mod3b <- update(mod3a, .~. + factor(lh):factor(func):factor(p))
anova(mod3b)
anova(mod3a)
Anova(mod3a)
 
# final model with one two-factor interaction
mod3c <- lm(pN^2 ~ factor(block)+factor(phy) + factor(lh)*factor(n)+factor(func)+factor(p), data=nitro[nitro$tag!=444,])
anova(mod3c)
 
# mixed effects model with crossed random effects: block and phy
mod3.lmer <- lmer(pN^2 ~ factor(lh)*factor(n)+factor(func)+factor(p)+(1|block)+(1|phy), data=nitro[nitro$tag!=444,])
anova(mod3.lmer)